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A diffusion-reaction model for the growth of bacterial 
colonies is presented. The often observed cooperative behav- 
ior developed by bacteria which increases their motility in 
adverse growth conditions is here introduced as a nonlinear 
diffusion term. The presence of this mechanism depends on a 
response which can present hysteresis. By changing only the 
concentrations of agar and initial nutrient, numerical integra- 
tion of the proposed model reproduces the different patterns 
shown by Bacillus subtilis OG-01. 

PACS: 87.10.-f-e, 87.17.Aa, 47.54.+r 



I. INTRODUCTION 

Some kinds of bacterial colonies present interesting 
structures during their growth [^-11|. Depending on 
the bacterial species and the culture conditions, colonies 
can exhibit a great diversity of forms. In general, the 
complexity of the growth pattern increases as the en- 
vironmental conditions become less favorable. Bacteria 
respond to adverse growth conditions by developing so- 
phisticated strategies and higher micro-level organization 
in order to cooperate more efficiently. Examples of these 
strategies are: the differentiation into longer-motile bac- 
teria, the production of extracellular wetting fluid, the 
secretion of surfactants which change the surface tension 
or the chemotactic response to chemical agents produced 
by bacteria . The experiments are usually made in a 
petri dish, which contains a solution of nutrient and agar. 
A drop of bacterial solution is then inoculated in the 
center of the dish. The growth conditions are controlled 
by the initial concentration of the medium components. 
The agar concentration determines the consistency of the 
medium, which becomes harder as the amount of agar in- 
creases, and the nutrient concentration controls the bac- 
terial reproduction. Depending on these two factors, the 
colony grows at a higher or lower rate, developing differ- 
ent kinds of patterns. 

In particular, colonies of the bacterium Bacillus subtilis 
OG-01 present a rich variety of structures ]^-[ll[|. Fig. 
|l| shows the morphological diagram obtained by Mat- 
sushita and co-workers [|lO|. They classified the colony 
patterns into five types, from A to E, whose main fea- 
tures can be summarized as follows. If the medium is 
very hard, i.e., with a high concentration of agar, bac- 
teria can hardly move and the colony essentially grows 



due to the consumption of nutrient and subsequent re- 
production. If the level of nutrient is also low (region 
A), the growth is controlled by the diffusion of the nu- 
trient up to the bacteria placed at the interface. The 
colony develops a ramified structure very similar to the 
patterns obtained with the diffusion-limited aggregation 
model (DLA) It takes approximately one month 

to cover the dish. If the initial agar concentration re- 
mains high and the nutrient concentration is increased, 
the growth is faster than in region A. The branches grow 
thicker until they fuse into a dense disk with rough inter- 
face (region B), similar to the patterns obtained with an 
Eden model ||l^. This structure needs 5-7 days to cover 
the disk. When the level of agar is decreased, which pro- 
duces a medium a little softer than in region B, and the 
level of nutrient remains high, the colony forms concen- 
tric rings (region C). This region is characterized by pe- 
riodic dynamics: for 2-3 hours the colony expands while 
the bacteria move actively (" migration phase" ) and then 
they almost stop for 2-3 hours ("consolidation phase"), 
during which the colony does not grow appreciably and 
the bacterial density increases due to reproduction. The 
crossover between the two phases is sharp. The periodic 
cycles of subsequent migration and consolidation phases 
create the pattern of concentric rings Accurate 
measurements show that in the growth phase there is a 
high concentration of longer and more motile bacteria, 
as a consequence of a differentiation process §. When 
a high level of nutrient is maintained, and the agar con- 
centration is decreased further, the colony spreads over 
the agar plate, and after less than 8 hours homogeneous 
disk of low bacterial density is formed covering all the 
dish (region D). In this thin surface, bacteria are always 
short and can move easily by swimming. By decreas- 
ing the nutrient concentrations for a semi-solid medium, 
the colony develops a densely branched pattern (region 
E) similar to the dense branching morphology (DBM) 
found in other systems |Q. The ratio of the width of 
the branches to the gap between them is constant over 
the whole colony. The colony grows quite fast, showing 
its main activity at the tips of the fingers, and covering 
the dish in less than 24 hours. The dynamics is related 
to both the consumption of nutrients and the bacterial 
motility. In general, when environmental conditions are 
adverse (low nutrient or hard surface), a higher level of 
cooperation is observed. 

The existence of a cooperative behavior seems to be 
determinant in the formation of the rings patterns of re- 
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gion C. The same kind of concentric rings has been found 
in experiments with other bacterial species. In the case 
of the bacterium Proteus mirabilis, the migration phases 
clearly involve the movement of differentiated swarmer 
bacteria (elongated and hyperflagellated) Similar 
ring patterns have also been observed in other non-living 
systems, like the Liesegang rings produced by precipita- 
tion in the wake of a moving reaction front [p_5l , or some 
experiments of interfacial electrodeposition |16[ | . In the 
case of the Liesegang patterns, it is well known that the 
distance between rings increases as t^^^, whereas in bac- 
terial and electrodeposition it is constant. 
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FIG. 1. Morphological diagram of patterns observed in 
colonies of Bacilus subtilis OG-01 as a function of nutrient 
concentration {C„) and the hardness of agar surface {1/Ca 
being Ca the agar concentration). Experiments were per- 
formed in a petri dish with a diameter of 88 mm. Taken with 
permission from 

Several models have been proposed to explain the vari- 
ety of patterns exhibited by Bacillus subtilis, as shown in 
Fig. |l|. DLA-hke patterns (region A in Fig. 

|l|) have been interpreted |^ as growth controlled by the 
diffusion of nutrients in the context of the DLA model. 
Ben- Jacob and co-workers proposed a communicating 
walkers model to describe some of the morphologies . 
This model reproduces the crossover between regions A 
and B by coupling random walkers to fields representing 
the nutrients. DBM-like patterns are also obtained by 
introducing a chcmotactic agent. Other kinds of mod- 
els are based on reaction-diffusion equations for bacterial 
density. The Fisher equation |17| can be used for repro- 
ducing the homogeneous circular morphology (region D) 
p8[ . Further developments were achieved by introduc- 
ing new elements to the Fisher model, such as a field for 
nutrient and nonlinear diffusion coefficients De- 
pending on the new elements introduced, these models 
reproduce some of the patterns of Fig. |l|. However, the 
ring patterns (region C) have so far eluded a satisfactory 



modelization. Although the model suggested in Ref. pc| ] 
can generate concentric ring patterns, they are rather 
different from those observed in experiments Q. In fact, 
dynamical cycles of consolidation and growth phases are 
not found. Finally, we must mention a model proposed 
by Esipov et al. |2l[| for the study of Proteus mirabilis 
colonies, which introduces a life-time for the differenti- 
ated swarmer bacteria. This model reproduces concentric 
ring patterns but does not explain why no periodicity is 
observed in other regions of the morphological diagram. 

In this paper, we propose a model consisting of two 
coupled diffusion-reaction equations for bacteria and nu- 
trient concentrations, where the bacterial diffusion coeffi- 
cient can adopt two different expressions, corresponding 
to two possible mechanisms of motion. The first is the 
usual random swimming performed by bacteria in liquid 
medium. The second is developed by bacteria in response 
to adverse growth conditions, and depends on their con- 
centration. Bacterial response is modeled as a global 
variable that can present hysteresis. Our model repro- 
duces the five morphologies observed in the experiments, 
including the ring patterns. 



II. MATHEMATICAL MODEL 

We consider a two-dimensional system containing bac- 
teria and nutrients. Both diffuse, while bacteria prolif- 
erate by feeding on nutrient. Let us denote by 6(r',T) 
the density of bacteria at time r and spatial position r', 
and by n(r',r) the concentration of nutrient. Then, b 
and n are in general governed by the following equations 
0@: 
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The function /(6, n) denotes the consumption term of 
nutrient by bacteria, and can be described by Michaelis- 
Menten kinetics 
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where k is the intrinsic consumption rate. For small n, 
the consumption rate is approximately linear in n and 
it saturates at the value fc/7' as n increases. Df, and 
Dn are the diffusion coefficients of bacteria and nutrient 
respectively. We assume that Z3„ is constant, but Di, can 
depend on nutrient and bacterium concentrations. 

As explained above, experiments show that in adverse 
conditions, bacteria can adapt themselves in order to im- 
prove their motility. In a soft medium and high nutrient 
concentration (region D of Fig. |l|), short bacteria can 
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swim randomly without difficulty, but in an adverse en- 
vironment (regions A,B,C and E) they need to develop 
mechanisms to become more motile. For intermediate 
conditions of semi-solid medium and sufficient nutrient, 
there are periods of fast growth (migration phase) and 
slow growth (consolidation phase) that lead to the con- 
centric ring patterns. 

The analysis of periodic rings suggests a dynamical 
scheme with hysteresis that can be outlined in the fol- 
lowing way: during the consolidation phase, the pop- 
ulation of longer-motile bacteria increases in order to 
overcome the opposition to the movement. When this 
population exceeds a certain value, enhanced movement 
becomes possible and a migration phase begins. Then, 
however, a progressive decrease in long-bacterial popu- 
lation ensues, until it reaches a minimum at which the 
" enhanced- movement mechanism" does not work. Then 
a new consolidation phase begins. Within this scheme, 
region D corresponds to a case where the maximum value 
is never reached (and therefore bacteria always move by 
usual diffusion) , whereas in regions A, B and E long- 
bacterial population does not fall below the minimum 
(and therefore always moves by the enhanced-movement 
mechanism). All these ideas can be introduced in our 
model by means of two basic points: 

(a) The diffusion coefficient Db{n,b) can take two dif- 
ferent expressions depending on the long-bacterial popu- 
lation. 

(b) The net production of long bacteria depends on the 
environmental conditions and also on the colony phase of 
growth. 

According to these ideas, we propose the following 
function for Di,: 



Db = D{di + d2b)fi, 



(3) 



where D depends on the concentration of agar, which 
is lower for harder medium. To take into account the 
inhomogeneities of the medium, we introduce a quenched 
disorder in D, which is written as D = Do(l+f(r')); ^{^') 
being a random term defined on a square lattice. From 
now on, Dq will be referred to as the diffusion parameter. 

The first term of Eq. (||) describes the usual diffusion of 
bacteria in a liquid medium. The second describes the co- 
operative enhanced-movement mechanism promoted by 
long bacteria. This second mechanism can be modeled 
by a diffusion coefficient that depends on the bacterial 
concentration. We multiply both terms by nutrient con- 
centration to take into account the fact that bacteria are 
inactive in the region where nutrient has been depleted. 
This dependence on n would not have been necessary if 
we had considered a "death" term in the equation for 6. 
The coefficients di and ^2 can adopt two different values 
(one of them zero) depending on the concentration of the 
long bacteria, as will be specifed below. 

Equations (|l|) with Eqs. (§)-(§) can be written in a 
simpler form as 
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di , d2 = d2 



At this point, we need to specify how to choose c?i and 
d2 depending on the population of long bacteria. In order 
to do this, we introduce a global phenomenological quan- 
tity W{t) that measures the amount of long bacteria. The 
evolution of this quantity should have a "creation term" 
that represents the transformation of short bacteria into 
long ones, and an "annihilation term" that represents the 
opposite transformation (septation) . It seems reasonable 
to assume that the creation term is directly dependent 
on the mean bacterial concentration, and inversely de- 
pendent on the level of nutrient (no) and on the diffusion 
parameter (adverse conditions, i.e. Dq and hq small, 
means a faster differentiation process). With regard to 
the annihilation term, it can adopt two possible values 
depending on the growth phase. The simplest equation 
that includes all these considerations can be written as 
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where A is a constant and c,- can have two different values 



> Cs). The quantity B, defined a,s B — 



'^b^ / ^b, is a measure of the mean concentration of 
bacteria inside the colony. We introduce the hysteresis 
previously pointed out by assuming that there are two 
limit values Wmax and Wmin for which: 

when W > Wmax then Ci = Cg , d2 — D2 , di — 0, 

(7) 

when W < Wmin then a = Cs , d2 — , di — 1. 

With a suitable choice of parameters A, Cs and Cg, and 
by changing only uq and Dq, we can obtain colonies that 
always move with one of the two types of diffusion, or 
colonies that periodically change from one type to the 
other. This will occur if Cs < XB /{hoDq) < Cg, which 
will give rise to the ring patterns. Although the bacterial 
response has been expressed in terms of the population 
of long-bacteria, other possible kinds of responses admit 
identical modelization. In this sense, our model is quite 
general. 
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III. NUMERICAL RESULTS 

We have numerically integrated Eqs. (^) with (^)-(0) 
in a square lattice of lateral size L = 600 using a 4-th 
order Runge-Kutta's method with mesh-size Ax — 0.5 
and time step At — 0.005. The system was initially 
prepared by assigning to each point a nutrient concen- 
tration n(r, 0) = no + r]{r, 0), rj being a uniform random 
number in the interval (—0.1,0.1), and a bacterial con- 
centration 6(r, 0) = 0, except in a small central square 
where 6(r, 0) — bo- The random term of the diffusion, 
^(r), takes a different and uncorrelated value in each box 
of side AAx. The random values are assumed to be uni- 
formly distributed in the interval (— e,e). The box size 
and the intensity e do not essentially affect the results. 

In all our simulations, we used the parameters bo = 0.7, 
7 = 0.5, £12 = 30, e = 0.4, A = 0.18, Cg = 2, = 1.6, 
Wmax = 3 and Wmin = 2. We reproduce the different 
morphologies observed in experiments by changing the 
values of the initial concentration of nutrients no and the 
softness of the media, related to Do- 
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FIG. 2. Bacterial colonies (top) and nutrient patterns 
(bottom) for a fixed value of Do — 0.01 and two values of 
initial nutrient: no = l(cs) and no — 5{b). They correspond 
to times t = 2500 and 75 respectively. 

In Figs. we present the results obtained for Do — 
0.005 and different values of rip. By increasing no we re- 
produce the crossover between regions A and B of Fig. 
|l|, from DLA-hke patterns (Fig. ||(a)) to a dense rough 
structure similar to that found with an Eden model (Fig. 



^(b)). All of them correspond to a situation in which, 
due to the small value of Dq, the creation term of Eq. (||) 
is greater than Cg, except at the very beginning. The re- 
sponse W can never decrease below the value Wmin , and 
therefore the colony will always grow with the enhanced- 
movement mechanism. In spite of this cooperative mech- 
anism, and because of the hardness of the medium, the 
effective diffusion DoD2bn is still small. The growth is 
mostly due to reproduction by feeding on the nutrient. 
For low level of nutrient, i.e- small ng, the colony growth 
is limited by the diffusion of these nutrients. It devel- 
ops branches, which are thicker as no increases. The 
prototype model that reproduces this kind of structure 
is the diffusion-limited aggregation (DLA) which is 
known to form a fractal pattern with a fractal dimension 
of dp = 1.71 Jl|,p2t. Experiments performed by Mat- 
sushita et al. in region A of Fig.l also show a fractal 
growth with dimension dp — 1.73 joj. We have ana- 
lyzed the fractal nature of the patterns obtained with 
our model, for Do = 0.005 and several values of the ini- 
tial nutrient, from rip = 1 (DLA-like) to rip = 5 (rough 
structure). We have calculated their fractal dimensions 
by using the box-counting method ||l|,^. In Fig. || we 
show, in a log-log plot, the number TV of boxes of size e 
that contains any part of the pattern, versus the size of 
the boxes. The slopes of the lines represent the fractal di- 
mensions. We observe that the cases that correspond to 
low nutrient have a fractal dimension of about dp ~ 1.73, 
showing good agreement with experiments. On the other 
hand, there is an abrupt change between these patterns 
and those that are not fractal {dp — 1). These last cases 
can be analyzed in terms of the roughness of their inter- 
faces. 
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FIG. 3. Fractal dimensions de- 

termined by the box-counting method, for Do = 0.01 and 
no = 1(0),2(0),3(D) and 5(A). Two lines of slopes 1 and 
1.73 are also plotted for comparison with experiments. 
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It is well known that Eden structures are not them- 
selves fractal, but their surfaces exhibit a self-affine scal- 
ing 23|. This implies that, for a long enough time, the 
width of the rough interface a scales with an exponent 
a as a function of the length of the interface I [a ^ 1°"). 
The roughness exponent for the Eden model is a = 0.5. 
Vicsek et al. [|3| analyzed experimental data correspond- 
ing to the region B of Fig. y. They concluded that these 
colony surfaces are self-affine with a roughness exponent 
a = 0.74. We have checked this point for our dense rough 
pattern (Fig. ^(b)) by measuring the width a for intervals 
of interface of length I. The results, as a function of ^, are 
presented in Fig. ^. In order to avoid additional effects 
derived from the radial growth of the colony, we have 
also performed a complementary simulation for the same 
parameters as Fig. ^b) but with a strip geometry. To do 
this, we have used a rectangular lattice of horizontal lat- 
eral size Lx = 600, with periodic boundary conditions in 
X direction, and taken as an initial condition for bacteria 
a horizontal line of length L^- The results for this case 
are also plotted in Fig. |^. For both circular and strip 
cases, we observe analogous behavior to that observed in 
experiments ]l3t . Our results show a linear region with a 
slope compatible with the experimental value a = 0.74. 




FIG. 4. Width of the rough interface as a function of the 
length of the interval in which it is measured. They corre- 
spond to the pattern of Fig. |(b) (□) and to a case with the 
same parameters but with strip geometry (A). A line of slope 
0.74 is also plotted for comparison with experimental results. 

With the aim of reproducing other morphologies of Fig. 
^, we now keep the initial nutrient fixed at the value 
no = 1 and increase the diffusion parameter Dq. Results 
are shown in Fig. |^(a)-(b). We observe a crossover from 
the DLA-like structure (Fig. ^(a)) to a dense branching 
morphology analogous to that represented in region E of 
Fig-i 



In Fig. 11(c)- (d), we present two snapshots obtained for 
a fixed value of the initial nutrient no = 5. They show 
how different kinds of patterns are obtained when Dq is 
increased: from the dense rough structure (Fig. ||(b)), to 
concentric rings (Fig. |^(c)) and homogeneous disk (Fig. 
^(d)). They correspond to the regions B,C and D re- 
spectively. Homogeneous disks are obtained when Dq 
and no are so high that the creation term of Eq. (^ is al- 
ways smaller than Cg. This means that the value Wmax, 
above which the enhanced-movement mechanism begins, 
is never reached, and bacteria move with the usual diffu- 
sion coefficient D^n. 




(c) (d) 

FIG. 5. Patterns obtained for no = 1 with Do = 0.05(a) 
and 0.1(6) and for no = 5 with Do = 0.1(c) and 2(d). They 
correspond to times t = 500, 300, 50 and 25 respectively. 

Ring patterns correspond to a narrow region of param- 
eters Dq and no for which the creation term of Eq. (||) 
takes a value between Cs and Cg . As explained in Section 
II, this leads to dynamics in which bacteria move alter- 
natively by usual diffusion Dqti (consolidation phase) or 
by the enhanced-movement mechanism DQD2bn (migra- 
tion phase). The two phases are clearly manifested in 
Fig. ^(a), where we represent the radius of the colony 
as a function of time. The pattern of concentric rings is 
a consequence of this dynamic behavior. In Fig.^(b) we 
plot the radial density profile, circularly-averaged, cor- 
responding to the ring pattern shown in Fig. ^(c). The 
maxima are formed in the positions where a consolidation 
phase began. To illustrate this point, we have pointed out 
in Fig. ^ the positions corresponding to the colony radius 
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at the beginning of each consohdation phase. 




FIG. 6. (a) Time evolution of the colony radius, for pa- 
rameters Do = 0.1 and no = 5; (6) Radial density profile 
corresponding to time t = 50 (pattern of Fig. ^c)). The 
positions where each consolidation phase started are pointed 
out. 

Numerically, our model also reproduces the ex- 
perimentally observed robustness of the growth-plus- 
consolidation period, which is barely dependent on 
changes in either nutrient or agar concentrations over a 
wide range. For high enough rtg, the value of the global 
quantity B approaches uq. In this limit, as can be derived 
from Eq. (p|) the period is given by 



T = Do{Wmax - Wmin) 
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which does not depend on no. Moreover, as a function 
of Dq, the period also maintains a rather constant value 
within a certain range (determined by parameters A, Cg 
and Cs) to increase sharply in the boundaries of the ring 
patterns region {Dq — > A/c^, — > ^/cg)- For equal 
period, the width of the rings increases with no. 



IV. CONCLUSIONS 

We have proposed a reaction-diffusion model for the 
study of bacterial colony growth on agar plates, which 
consists of two coupled equations for nutrient and bacte- 
rial concentrations. The most important feature, which 
introduces differences from previous models, is the fact 
that here we consider two mechanisms for the bacterial 
movement: the random swimming in a liquid medium, 
and a cooperative enhanced movement developed by bac- 
teria when the growth conditions are adverse. The two 



mechanisms are introduced in our model by means of 
a diffusion term with two different expressions which de- 
pend on the bacterial response to the environmental con- 
ditions. This response is modeled as a global variable 
that presents hysteresis depending on the conditions of 
the medium. The inhomogeneities of the agar plate have 
been taken into account as a quenched disorder in the 
diffusion parameter. 

We have shown that, simply by changing the parame- 
ters related to the hardness of the medium and the ini- 
tial nutrient, our model reproduces all the patterns ob- 
tained experimentally with the bacterium Bacillus sub- 
tilis: DLA-like, dense-rough disk, DBM-like, ring pat- 
terns and homogeneous disk. We have calculated the 
fractal dimension of the DLA-like structures and the 
roughness exponent of the rough disk surface, obtain- 
ing results in good agreement with experiments. The 
ring patterns have been obtained for intermediate values 
of agar and high nutrient. In this region, the bacterial 
response presents hysteresis and the two mechanisms of 
motion work alternatively, leading to cycles of migration 
and consolidation phases. The duration of these cycles 
is roughly constant for different values of nutrient and 
agar concentration over a wide range. This periodical 
dynamics generates patterns of concentric rings. 

In summary, the model proposed satisfactorily repro- 
duces the whole experimental morphological diagram. It 
represents a first attempt at describing the response of 
bacteria to adverse growth conditions and, in certain con- 
ditions, their ability to improve their motility. Further 
refinements could be made. The bacterial response, here 
described as a global variable W{t), could be considered 
in a more realistic way by introducing a coupling term 
in a local version of Eq. (o) for a field W{r, t). However, 



preliminar studies with such a model 
essential features previously describee 



shows the same 
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